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SUMMARY 


The Inverse method, with the shock wave prescribed to be an 
elliptic cone at a finite angle of incidence, is applied to calculate 
numerically the supersonic perfect -gas flow past conical bodies not hav- 
ing axial symmetry. Two formulations of the problem are employed, one 
using a pair of stream functions and the other involving entropy and 
components of velocity. A number of solutions are presented, illustrat- 
ing the numerical methods employed, and showing the effects of moderate 
variation of the initial parameters . 


INTRODUCTION 


The first three-dimensional high-speed flow problem to be solved 
exactly was the flow past a circular cone. Pioneers in this work were 
Busemann, and Taylor and Macoll. (See refs. 1 to 3*) Reference 4 
gives numerical solutions for circular cones using the theory developed 
in reference 2. A perturbation technique was devised by Stone, refer- 
ence 5, for circular cones at small angles of yaw. The method of Stone 
was applied in reference 6 to obtain approximate solutions for cones at 
small angles of yaw, and was extended in reference 7 for large angles 
of yaw. A somewhat similar numerical analysis, the linear characteristics 
method, was developed by Ferri for analyzing the flow past yawed circular 
cones (ref. 8) and slightly noncircular cones (refs. 9 a-nd 10). Tables, 
based on the method of Ferri, are available (ref. ll) for calculating 
such flows. Of the various analyses referred to above, the only exact 
numerical solutions, based on the full inviscid equations of motion, 
are the ones tabulated in reference 4. They were obtained by numerical 
integration of the equations of motion specialized for symmetrical 
conical flow. 

There are various ways that an exact (numerical) attack on a high- 
speed flow problem, employing the full inviscid equations of motion, 
might be organized. One is to specify conditions at the body and then 
to integrate the equations of motion outwardly toward the shock wave. 
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This is the manner in which the problem was posed for the integrations 

tabulated in reference 4. Another way is to prescribe conditions at 

the shock wave and then to integrate inwardly t award the body. This has 

been called the "inverse" or "marching inward" method. The capabilities 

of this method are exemplified in the work of Van Dyke and Gordon , 

reference 12, who presented a catalog of solutions of flows past blunt 

bodies for a large family of preassigned shock shapes. In the work of 

Radhakrishnan , reference 13, the marching inward technique was used to 

find the conical body behind a circular conical shock wave at a large 

angle of yaw. In reference l4, Briggs employed an analogous approach 

to compute bodies that support elliptic conical shocks. A 

The work of reference l4 is extended here zo include circular or 
elliptic conical shock waves at finite angles of attack or yaw. Two 
implementations of the inverse method are employed, and each has been 
programmed for an electronic computer. The first formulation, which 
was used in reference l4, is in terms of pressure, density, and a pair 
of stream functions , one of which vanishes at the body. (A discussion 
of the use of several stream functions in the analysis of three- 
dimensional flows will be found in reference 15 ) The second formulation 
of the problem is in terms of velocity components, pressure, density, 
and entropy. The body is located by calculating the position of the 
entropy layer, shown by Ferri (ref. 8) to exist on the surface of the 
body in the case of nonsymmetrical conical flow. The computations of 
reference 13 were based on this latter approach. (An interesting 
discussion relative to the thickness of the af orementioned entropy layer 
is given by Cheng, ref. lo.) 

The flow field behind an elliptic conical £ hock wave is analyzed 
in detail in the present report, both the stream-function and velocity- 
entropy formulations of the problem being used. The resulting body and 
surface pressure is given, along with plots of 1 he variation of the 
velocity components and entropy in the flow field. A diagram showing 
the stream lines in the crossflow is also presented. Comparisons are 
shown between results of wind-tunnel tests on elliptic cones and numerical 
calculations with initial data based on s chile rc n pictures of the shock 
waves in these tests. A comparison with the apjroximate theory of 
reference 6 is made using a machine result calculated with a circular 
conical shock wave at a finite angle of yaw. Finally, the results of a 
number of cases are presented wherein systematic variations of shock 
geometry, free-stream Mach number, and ratio of specific heats were made. 

These are given in the form of body shape and surface pressure. 


SYMBOLS 


a i,bi 


semimajor and semiminor axes of elliptic cross section of 
coordinate surface of constant 0 = 0 O (eqs. (4)) 
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semimajor and semiminor axes of elliptic cross section of 
coordinate surface of constant cp = cp Q (eqs. (^)) 

functions defined by equations ( 8 ) 

stream functions of a three -dimensional flow 

specializations of f and g for the present conical flow 
problem 

coefficients of the metric of the sphero-conal coordinate 
system 

constants in the sphero-conal coordinate transformation 
Mach number 

number of points taken in the numerical integration 

function defined by equation (4l) 

pressure referred to PoqU^ 2 

sphero-conal coordinates (see fig. 2) 

entropy, JL 
P? 

velocity components in the sphero-conal coordinates 

velocity components in Cartesian coordinates 

free -stream velocity 

velocity vector 

Cartesian coordinates 

ratio of specific heats 

coordinate surface that coincides with the shock wave 

increments of 6 and cp that define the mesh size in the 
numerical calculations 

density referred to free -stream value 



Subscripts 
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conditions on the body surface 
nth extrapolation 
normal component (eq. (23)) 
conditions at the shock wave 

differentiation with respect to variables r,0,9 
differentiation with respect to variables x,y, z 
free -stream conditions 

Superscripts 

(n) nth extrapolation 

r differentiation with respect to the variable G 

THE SPHERO-CONAL COORDINATE SYSTEM 


b 

n 

N 

0 

x,y,z 


The numerical procedure employed here requires that the shock 
wave be a coordinate surface. A system of ortnogonal coordinates 
r,Q,ty which meets this requirement is given in reference 17 . The 
transformation from Cartesian coordinates is 

x = r cos 0 Yl-k* 2 cos 2 cp 

y = r sin 0 sin 9 

z = r cos 9 n/ l-k 2 cos20 

k 2 + k 1 2 =1 } 0 < 9 < 2jt | 

where x,y,z are the Cartesian coordinates. The element of length in 
the sphero-conal coordinates system is defined by the relation 



ds 2 = h x 2 dr 2 + h 2 2 d0 2 + h 3 2 d9 2 


( 2 ) 
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where 


h n 2 = 1 


h 2 2 


= r 2 


k 2 s in 2 0 + k ,2 sin 2 cp 
l-k 2 cos 2 0 


hr, 2 = r 2 


k 2 sin 2 0 + k t2 sin 2 cp 
1-k* 3 cos 2 q> 


(3) 


The coordinate surfaces are: 
x 2 + y 2 + z 2 = r 2 


x 2 


yf 


bi 2 


a i e 


= 0 


for fixed r = r Q 
for fixed 9 = 0 o 




2 2 i 

— 2 + ~ z 2 = 0 for fixed cp = cp Q 

2 _ l-k 2 cos 2 0<-, \ (*0 

1 k 2 cos 2 0 o 

b^ = tan 2 0 o a x > bi 

2 _ l-k ,2 cos 2 cp 0 
&2 k ,2 cos 2 9o 

b 2 tan 2 cpQ a 2 = b 2 

J 

Portions of the conical surfaces that are represented by the second and 
third of equations (4) are shown in figure 1. The shock wave of the 
present problem will be taken as a constant 0 ( =0 O ) surface . 


The components u f v, and w of velocity in the sphero-conal 
coordinate system are in the directions of increasing r , 0 , and cp, 
respectively . Figure 2 shows curves of constant 0 and cp as projected 
onto the y-z plane- The velocities v-]_ and w x are in the directions 
of increasing y and z, respectively, and the velocity ui, not shown 
in figure 2, is in the direction of increasing x. 
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The relationship between Cartesian and sphero-conal components of 
velocity is given by the equations: 



THE EQUATIONS OF MCTION 

The equations of motion for flow of ar ideal compressible fluid 
are, in vector form, 

div(pV) = 0 (continuity) (6a) 

p(V.grad) V + grad p = 0 (momentum) (6b) 

V.grad(p/p7) = 0 (energy) (6c) 

The velocity V has been made dimens ionles s with respect to free -stream 
velocity Uqq, density p with respect to j'ree-stream density p^, and 

pressure with respect to the quantity PoJJ^ 2 - When transformed into 
the sphero-conal system, equations (6) become 

A M P v -)_ + B d( p y -), + p ( 2 u+ACv+BDw ) = 0 

de 6q? 


(7a) 
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Av 


Ay 


Av + ^ " (v 2+ “ 2 ) ■ o (7b) 

Se + Bw 5 ” ACv2 + uv * BD ™ + p te = 0 (7= ) 

§ + B “|- BD ^ + ““ + AC ™ + f| = o ( 7 a) 

Av ♦ BW Zgil - 0 (7e) 


Note that the assumption oi conical flow has "been made in the writing 
of equations (7) and that the following abbreviations have been employed: 


A 2 = 


l-k 2 cos 2 0 _ r 2 

k 2 sin 2 0 + k ,2 sin 2 c p h 2 2 


B 2 = 


C = 


l-k T cos 2 cp 


k 2 sin 2 0 + k t2 sin 2 cp h 3 2 


k 2 sin 0 cos 0 


k 2 s in 2 0 + k T 2 s in 2 cp 




k T2 sin cp cos 9 
k 2 sin 2 0 + k f2 sin 2 cp 


D = 


( 8 ) 


Three-dimensional flows such as those under discussion can be 
described by a method utilizing two functions which are a generalization 
of the familiar stream function of two-dimensional and axisymmetric 
flows. This technique was used in reference 1 ^ and a concise description 
of the use of stream functions in the analysis of three-dimensional flow 
problems can be found in reference 15. 

If two functions^ f and g, are chosen such that 


pV = (grad f ) X (grad g) 


(9) 



then equation (6a) is satisfied, and f and § represent stream ^surfaces 
whose intersections are streamlines . This can he demonstrated in the 
following manner. 

Equation (9) has the three Cartesian components 

pu = fyg z “ ^z§y 

P v = ^z&x “ f x&z 

pw = f x gy - fyg x 

where subscripts x,y,z indicate partial dii 
equations (10a) it follows that 

u:v:w = (fyg z - f zgy ) : (^zSx " 

Equation (10b) is a solution of the two equa^ 

V-grad f = 0 

V-grad g = 0 

Equations (ll) are valid if V, the velocity vector, is normal to grad f 
and grad g. Grad f and grad g are the normals to surfaces, f and g, 
and so V must lie in the surfaces f and g. It can be said, then, 
that f and g surfaces contain the streamlines, and that the inter- 
sections of these surfaces are the streamlines . Thus it is quite 
reasonable to think of f and g as stream functions. 

The energy equation, equation (6c), now takes the form 


= S(f,g) (12) 

which states that entropy is constant along streamlines . The momentum 
equation, equation (6b), can be written in "terms of V, p, and S; that is, 

p(V-grad)\f + grad[p?S(f ,g)] = 0 (13) 

The velocity components u,v,w can be found in terms of the stream 
functions f and g by expanding equation ( c >)* Th e result, in the 
sphero-conal coordinate system, is 


(10a) 

A 

3 
8 

f erentiation . From s 


):(f X gy - %g x ) ( 10b ) 

ions 


( 11 ) 
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A 


U pr 2 ^ f 0 S( P " f 9 s 0 ^ 

g 

~ pr ( ^cpSp “ f r gcp ) 

A 

w = ^r ( f rS0 - f 0gr) 


) (1*0 


Here f and g are functions of r,d,cp and the subscripts indicate 
partial differentiation. Equations (l4) are general for three- 
dimensional flows and are not limited to conical flows. 

• — ► 

For conical flow V, p, and S are not functions of r. Thus it 
is clear from equations (l4) that the product fg is homogeneous of 
order 2 in r. In order that f will reduce to Stokes 1 stream function 
in the axisymmetric case and that the tody be described by f = 0, it 
is convenient to make the definitions 

f(r, 0 ,q>) = r2F(e,cp) 

( 15 ) 

g(r,9,cp) = G(0,q>) 

Now S cannot be independent of r unless it is also independent of 
f . It can be concluded, then that 

S(f,g) = S(G) ( 16 ) 


for conical flow. 


Equations (l4) may now be formally written in terms of F and G; 
that is , 


u = 


AB 


(F^Gcp - FcpGg ) 


^ 2B ™ 
v p FG P> 


(17) 


2A 

P 


V = ^ FG 0 


The components of 


grad[p 7 S(G) ] are 
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> 


(18) 


in the r,0,cp directions. The equations of motion in terms of F and 
G are: 


0 momentum: 


[7P 7+1 S - (2BFG 9 ) 2 ] ^ 


= (2BF) 2 (G 0 G cpcp -G C pG0 C p) - (p 7 + 1 S) (-g-J G q + 2B 2 FGcp(Fc p G 0 -F 0 G cp ) 


+ (2F) 2 G 0 (d-B 2 -f) G 9 + (2F) 2 (B 2 G cp 2 + A 2 Gq 2 ) 


<p 


(19) 


cp momentum: 


(2AF) 2 G cp G 00 = 2A 2 EG 0 (F cp G 0 -F 0 G cp ) + (2APG 0 )‘ 


J 9 cp " 


- G fl ( D + ^ + G m — 


p 


(2F) 2 G 9 (CG 0 + B 2 DG(p) + p? +1 S ( 7 ^ G, 


J cp 


p cp . s’ 


p T S 


( 20 ) 


r momentum: 


A 

3 

8 

9 


(abg 9 ) 2 f 00 


= (AB) ; 


F0G 0 Grp<p - (F 0 G C p+F C pG 0 )G 0C p + F C pG C pG 00 - G 0 2 Fq X p + 2G 0 G 9 F 0C p 


- 2F(A 2 G 0 2 +B 2 G cp 2 ) - (F C pG 0 -F 0 G C p) 


3 2A C-C+A e 


p e 


-'cp 


- A 2 ( 2B 2 D-D+B 2 j G 0 


( 21 ) 
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It has been stipulated that F is to vanish on the surface of the 
body. It would appear, then (see eqs. (lj)), that the velocities v 
and w must also be zero there. This is not generally the case, so 
the derivatives G<p and Gq must become infinite on the body. The 
singular nature of the derivatives of G at the body is evidence of 
the previously mentioned singular behavior of the entropy near the body, 
since S and G are functionally related (see eq. (l6)). 


THE INITIAL CONDITIONS 


The shock wave is taken to be an elliptic cone of the family 
described by the second of equations ( A ) . The value of 0 at the shock 
cone is designated 0 O . The equation of this conical surface is 


r U,y,z) = x 2 - - .Jn. = 0 

h i a i 2 


( 22 ) 


The component of the free -stream Mach number normal to the shock wave is 

+ «z( 


% “ 


(- 

"\Sx 


MsdS-J + 


Sr 

dy. 


Sr 

Sz 


Sr 

Sx 


Sp 

.Sy. 


Sr 

Sz 


( 23 ) 


where M x ,My-,M z are the x ; y, and z components of the free -stream 
Mach number. When rewritten in terms of the sphero-conal coordinates, 
equation (23) becomes 


^ = sin 0 o cos a cos p / (l-k 2 cos 2 9 9 Xl-k' 2 cos 2 cp) 

J k 2 sin 2 0 o + k^sin 2 ^ 


-sin 0 o sin cp cos a sin 3 




l-k 2 cos 2 0 o 


k 2 sin 2 0 o + k ,2 sin 2 cp 


- cos cp 


sin a 


k 2 sin 0 o cos 0 O 
N/k 2 sin 2 0 o + k !2 sin 2 cp 


(24) 
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The angles of attack and yaw, a and |3, and the x, y, and z components 
of Mach number are shown in figure 3- Note that the axes are fixed in 
the shock wave and that their origin is at the cone vertex. 

Expressions are required for evaluating velocity components, density, 
pressure, and the stream functions just behind the shock wave. The 
oblique shock relations (see, e.g., ref. l8) along with equation (24) 
and equations ( 17 ) are utilized in the writing of these expressions. In 
the following paragraphs the equations that apply only at the shock will 
be clearly indicated so as to distinguish them from such general relations 
as may appear in the analysis. 

The velocity components u and w lie in a plane tangent to the 
shock, just ahead of the shock. They therefore undergo no change through 
the shock, so that at 9 = 0 O 


u 


cos 0 o cos a cos 


p n/ 1-k' 2 cos 2 cp 


5 


+ 


sin 0 o sin cp cos a sin (3 + cos p sin a 




-k 2 cos 2 0 o 


(25) 


w = cos 0 o cos a cos 


k' 2 


sin. cp ccs cp 


Vk 2 sin 2 e o + k^sin^ 


- sin 9 sin cp cos a sin 



1-k' 2 cos 2 cp 
k 2 sin 2 0 o + k' 2 sin 2 cp 


- sin cp 


sin a 


(l-k 2 cos 2 0 o ) (l-k 1 2 cos 2 cp) 
k 2 sin 2 0 o + k' 2 sin 2 cp 


( 26 ) 


The velocity v and the density and pressure are found from the shock 
relations to he 


_ (y-i)Mn 2 + 2 

(7+l)MjjMoo 

(y+l)M M 2 
(y-l) M N 2 + 2 


(27) * 

( 28 ) 


00 Co > 
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27 % 2 - (7-1) 
7( 7+1 )Moo 2 


( 29 ) 


where M-^ is represented by equation (24). 

If equations ( 27 ) j ( 28 )., and the first of equations ( 17 ) are 
combined > the result is 


% 

PV ' ~ Hx, ~ 


(30) 


If the second equality is taken in equation ( 30 ), and if equation (24) is 
inserted for Mjj an expression is obtained which can be solved for 
G<p, valid, of course, only at 0 = 0 O ; that is. 


JL_ 

2F 


Gcp _ pv ( cos ©ocos a cos p *] l-k 2 cos 2 6 


- cos 0 osin cp cos a sin (3 


I l-k 2 cos % 0 
1 -k* 2 cos3cp 


- cos cp sin a 


k 2 sin 0 o cos 0 O \, 
•J 1-k' 2 cos 2 cp ^ 


( 31 ) 


at 0 = 0 O . The form of the functions F and G at the shock is still, 
to some extent, arbitrary. If it is asked that G be equal to cp on 
the shock, then this can be accomplished by setting 


F = — ( sin 0 o cos a cos p J l-k 2 cos^9 c 


- cos 0 o sin cp cos a sin 


in p l ~ — 

J 1 -*’ 


-k ^cos^p 
2 cos 2 cp 


- cos cp sin a 


k 2 sin 0 o cos 0 O 1 
7l-k ,2 cos 2 © ^ 


(32) 



The streamlines , upon 'which entropy is constant, converge near the "body 
and come together at B. Thus all values of entropy exist at B, a 
singular point. The point E is also singular hy reasons of symmetry. 
The entropy therefore undergoes rapid change in the vicinity of the "body, 
thus leading to the notion of an entropy layer. The streamline emanating 
from the shock at C impinges upon the body an A, and follows the body 
surface around to B. The velocity v is normal to the body at points 
A and B, and so it is zero there. The body can be completely determined, 
then, by plotting the points on AC and BD where v vanishes and by 
noting the points where curves S versus cp (0 = const.) cross the line 
S(0 = 0 O , cp = 0) = const. 


THE NUMERICAL ANALYSIS 


The conical shock problem is solved here by two methods, one using 
the stream-function formulation of the equations of motion and the other 
using the velocity-component formulation. The two procedures are pro- 
grammed for electronic computers, and for convenience they are termed the 
"stream-function program" and the "velocity-entropy program." Descriptions 
of the computer programs are given in the following sections, along with 
a brief discussion of their specific application and limitations. 


The Stream-Function Program 


jP]^xs program has been set up for an Electro— Data 205 Electronic Data 
Processing Machine. The program is in most ways identical to the one with 
which the calculations were performed in reference l4, with the extenstion 
here to include angles of attack or yaw. Equations (19) to (2l) , along 
with equations (25) to (28) , (33) "to (36) , anc. (38) to (42) , define the 
problem. 

For a given case, values are prescribed for 


bi/ai 

M» 

7 

a 


shock-wave geometry 

free-stream Mach number 
adiabatic exponent 
angle of attack 


angle of yaw 
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AO increment in d } a positive quantity 

n number of points 


The increment Acp is found from the tabulation below, which is 
incorporated in the computer program. 


a p Acp 


Range of cp 


0 0 rt/ ^n 

0 /0 Ttj n 

/0 0 jt/n 


(Acp/ 2) < cp < (rt/2) - (Acp/ 2) 

-(rt/2) + (Acp/ 2) < cp < (at/ 2) - (Acp/ 2) 
(Acp/ 2) < cp < 3t - (Acp/ 2) 


In practice 15 or 20 points are used, but as few as 10 and as many as 60 
have been used in some cases. The size of A 6 is chosen so that 8 or 
10 increments are taken to get to the body. Cases have been run with as 
few as 2 and up to 50 or 60 steps to the body. No particular difficulty 
is encountered in using the very small increments. This is in contrast 
with the calculations of reference 12, where an inherent instability limits 
the number of steps to the body. 


The calculations are arranged according to the outline below. 


Step la .- Compute initial values of p, F, G, F^ , and 
at the shock wave using the input values of 0 o , bi/ai, etc. 

Step l.» Compute the derivatives p , V ^Ccp^ %cp ' 
F^, and G^ numerically. Compute S and S ! /S. Solve the 
equations of motion for Pq, Fqq , and G^. Calculate the 
"externally iterated" value of p, designated P±-t> using the 
formula 


P 


(n) 

it 



1) 


Ae f (n) (n - 1) 

T p e p e 


(At the shock, p^+ = p.) It is termed "externally iterated" 
in that is is computed at each point using only the informa- 
tion at hand, and it is not used in the subsequent extrapo- 
lations. The pressure p is computed with the relation 


p = p-j^S. Pressure, density, and other pertinent flow-field 
data are now read out of the computer. 


(43) 


Step 2 .- Extrapolate p, F q, and Gq to the coordinate 8 n+1 
using the formulas 
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p(n + i) = p (n) _ A0p( n ) 

0 

4 " + i) = F (n) _ ^(n) , 

4 n * 11 - 4 n) - ^44 

> 

Take these values of Fg + and G^ 11 + along with the 

previously calculated results and G^ and extrapolate 

F and G with the equations 


i<(n + i) _ 


- Tf^ 11 ^ 4- R 1 ' + 


,(n + l) 


££ G^ n) + G^ r ' + ^ 
2 9 e 


Step 3 .- If any values of the stream function F have 
become negative at 9 = 0 n + x go to step 4. Otherwise 
repeat steps 1, 2 y and 3* 

Step 4 .- Some values of F are negative at e = 0 + 1 . 

Recall the flow data for 0 = 0 n and extrapolate to the body 
with the following equations. The subscript b indicates 
results on the body. 


A0 b = Vr 

f 1 ^ 


0b “ ^n “ ^b 


sin cp sin 0^ 


/ 2 2 
cos Q^sj 1 - k * cos 0 


cos cp 

COS 0-j^M 


fl - k 2 cos 2 0i ) 
1 - k T cos cp 
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P b = P (n) - Aa bP ^ n) 





+ 



(n) (n) i 

G e 


m 


Note that y b and z b are actually ratios (y/x)b and ( z/x)b, 
and are , as such, projections of the body coordinates r, 
t*b> 9 onto a plane x = 1. These coordinates, along 

with ph , are now read out of the computer and the case is 
finished. 


The order of procedure discussed in the foregoing paragraphs can be 
seen graphically in figure 5, which is a simplified flow chart of the 
computing program. 

Note in step 3 of the foregoing outline that the integration is 
stopped as soon as the stream function F changes sign anywhere on the 
coordinate ellipse 0 n + The body is then located by linear extrapo- 
lation of the data at to the place where F vanishes. In cases 

where the ellipse & n is markedly different in shape from the body, the 
extrapolated distance, A0 b , may be excessively large, resulting in an 
inaccurate representation of the body. This situation occurs for very 
large values of cc or $, or for small values of b^^a^, or for small Mach 
numbers. The computer program was not readily amenable to modification 
to overcome this difficulty, which in general prevented close analysis of 
the entropy layer on the body. A number of cases were calculated using 
this program, and limits on parameter values were generally established 
from comparison of results with calculations using the velocity— entropy 
program. The latter program has been set up so as to overcome the 
difficulty discussed above, and it is described in the following section. 


The Velocity- Entropy Program 


This program was set up for an IBM 704 electronic computer for the 
purposes of studying the entropy layer on conical bodies, and as an 
independent check on the stream- function program. The velocity formula- 
tion of the equations of motion, equations (7) , is used here. This 
formulation is numerically less involved than the stream- function formu- 
lation. Only first derivatives must be taken numerically in the velocity 
formulation, for example, whereas second derivatives are required in the 
other. The velocity-entropy program is the more versatile of the two in 
that the range of cp is arbitrary and need not be precisely it/2 or jt 
radians. It is the more difficult of the two to use, however, in that 
the body is found by an analysis of the entropy, which cannot easily be 
programmed for the computer. It is not intended as a production program 
and only a few cases have been analyzed in detail with its use. 
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The initial shock-wave data are found from equations ( 25 ) to ( 29 ) 
with prescribed values of Q 0 , bi/ai, 7 , a, (3, A G , and n. The number 
of points n may he in the range 7 < n < 50- (The lower limit 7 is 
dictated by the numerical differentiation scheme that is employed here.) 

The range of cp is also specified along with 6 Q , bx/ax, etc., by pre- 
scribing values of Lx and L 2 . These two quantities are angles and are 
put into the relation Acp = (L 2 - L 1 )/n. The values of cp are then 
taken at intervals Acp in the range (Lx + Acp/ 2) cp < (L 2 - Acp/2) . For 
a=[3 = 0 it is usual to take Lx = 0 and L 2 - 90°, but different values 
may be used if desired. An outline of the commuting procedure is presented 
next . 


Step la .- Compute initial values of u ; v ; ¥ ; p, and p at 
the shock wave using the given values of 0 Q , bx/ax, etc. Compute 
s = p/pn 

Step 1 .- Calculate numerically the derivatives ru , v^, w^, p^,, 
S^. Calculate ug, Vg , w g, pg , and Sq from the equations of 
motion. Compute p.^ using equation {43) . Compute pressure p 
using the formula p = Pit^S* Read current values, at 0 = 0 n , 
of u, v, w, p, p, and S out of the computer. 


Step 2 .- Extrapolate velocities 
coordinate using the formulas 

(n + 1) (n) 

u = u v 

v (n + = v (») 

w(n + ^ 

p (n + 1) _ p (n) 
g(n + 1) _ g (n) 


p, and S to the next G 


“V 

ASpta 

A0S^ n ‘ 




> 


J 


(^9) 
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Step 3 >- any values of the velocity v have become posi- 
tive at G = 0n+i £° 4o step 4. Otherwise -epeat steps 1, 2, 
and 3* 

Step 4 .- Some values of v ( n+1 ) have become positive. Recall 
the flow-field results for the station G = On and continue the 
integration using one of the two options A and B. The option 
choice is specified along with the initial data G a , bx/ax, etc. 


< 00 co 
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Option A.- Calculate a new value for A0 . The new A0 
is some fraction, say 1/5, of the original value. The integra- 
tion is continued now (steps 1 and 2 ) , and instead of using all 
n data points at once, several smaller segments are used. The 
segments may overlap, and the integration proceeds on each 
segment until all values of S are greater than the value cal- 
culated to exist on the "body, or until a pre-set number of 
steps have been taken, or until the calculation breaks down as 
evidenced by negative S or p values. The size of the small 
segments, the number of points between the first points of 
succeeding segments, and the size of the new A 0 are all 
specified with the initial parameters . 

Option B .- This process starts at 0 = 0 n in the flow 
field and uses a reduced value of A 0 , as in option A. It 
depends on the entropy variations being of the form shown in 
figure 6 , a situation that prevails for elliptic cone shocks 
at a = J3 = 0 in the region 0 < cp < 90°, or circular cone 
shocks at finite angles J3 in the range -90° < cp < 90°. 

Observe in figure 6 that the derivative is negative over 

most of the flow field, but that zero values begin to appear 
as the integration proceeds. The points of zero slope are 
labeled S± and S 2 • If the data to the right of point S 2 on 
the 0 n + 2 curve are used in extrapolating to the coordinate 
0 n + 3 ; numerical difficulties arise due to extreme values of 
Sq) in this region, and due also to the fact that the velocity 
v passes through zero near S = 0. Divisions by very small 
values of v can introduce large errors in the subsequent 
calculation of S, p, etc., from the equations of motion. These 
same difficulties appear in data to the left of points S x . The 
technique in option B is to continue the integration from 0 n , 
but to retain, at the extrapolation step, only the data lying in 
the region between points Si and S 2 on the entropy curve. The 
integration is continued in this manner until erroneous negative 
values of p or S occur within the range cp(Si) to cp( s 2 ) , or 
until all values of S are greater than the body value. 

The body is found with either option by converting the cp - 0 
coordinates of places where the entropy curves cross the line of 
constant entropy that represents the body entropy to y^ and 
using equations (47) • 

No flow chart is presented for the velocity-entropy program. It is 
identical in concept to the one shown in figure 4 for the stream-function 
program. 

As has been indicated, the numerical data near the body are subject 
to the ravages of division by zero and abrupt changes in slope of quanti- 
ties whose derivatives must be obtained numerically. These problems occur 
to the largest extent in the region nearest the entropy singularity, and 
have a physical basis in the fact that the entropy does indeed change very 



rapidly in a thin region near the body, and most rapidly near and at the 
singular point . There does not seem to be any general numerical procedure 
suitable for machine programming by which the calculations can be performed 
with assurance in this region. Results can usually be obtained over the 
first portion of a body, but only with difficulty over the region ending 

at the singularity, even with the devices of options A and B. In a given 

case it is usually necessary to recompute several times, using both options 
A and B, and with various lengths of segment size in option A in order to 

get a reasonably complete picture of the entropy variation in the region 

near the entropy singularity. 


PRESENTATION OF RESULTS 
Preliminary Calculations and Comparisons 


Preliminary calculations for elliptic ccnical shock waves at zero 
angles of incidence are given in reference 1+. The main purpose of these 
calculations was to study the stability and convergence of the numerical 
process, and the effect on the results of variations of 7 and It 

was found that the calculations are stable and convergent for small incre- 
ments A0 . The bodies that support elliptic conical shocks are themselves 
somewhat flatter than elliptic cones. The change in body shape and surface 
pressure with changing Mach number is as might have been anticipated. The 
limiting case of 7 = 1, M*, — °o corresponds to the so-called Newtonian 
plus centrifugal theory. A machine result calculated with 7 = 1.001 and 
Mm = 10,000 is compared with the simple Newtonian theory, which states that 
the surface pressure is proportional to the square of the sine of the local 
body angle. Laval, reference 19, has calculated pressures on conical 
bodies in the limit 7 = 1, Moo = <». He presents a result for the case of 
do - 30°, bi/ai = 0.6, and a = |3 = 0, the suae case used for the comparison 
in reference l4. These results are brought together in figure 7. 

Initial comparisons of calculations using the two computer programs 
described in the present report were made wish the case of 0 O = 30°, 
bi/ai = 0.6, = 6, 7 = 1.4, and a = p = 0. The computed body shape 

and surface pressure are shown in figure 8. Observe that the results of 
the two computer programs differ by a small amount. Another comparison 
of the two computer methods is shown in figure 9 for the case of a cir- 
cular shock wave, 0 O = 30°, at an angle of yaw of 20°, and with 7 = 1.4 
and Mao = 10. This is the case reported in reference 13. Again, there is 
a slight discrepancy between the results of the two machine procedures. 

The velocity-entropy calculation agrees closely with the result of refer- 
ence 13, except on the leeward side, which is remarkable when it is 
realized that Radhakrishnan' s work was carried out on a desk computer. 

The differences in the machine calculated results can probably be ascribed 
in part to differences in programming procedures, and to the numerical 
difficulty associated with the body extrapolation process of the stream- 
function program. The differences are not great, however, but the velocity- 
entropy program is presumed to give the more precise results. 
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The variations of entropy and velocity components with the coordinate 
cp for the case shown in figure 8 are presented in figures 10 to 13- The 
streamline pattern in the crossflow is shown in figure It. 

Stocker has calculated conical “bodies by a two- stream- function 
technique wherein the stream functions are independent rather than depend- 
ent variables, No report has been published by Stocker , as yet, but the 
results of one of his calculations, as given in reference 20, are shown 
for comparison in figure 8. His results agree well with calculations using 
the velocity-entropy program. 

Mauger, reference 21, has applied complex variable theory and the 
method of Garabedian (ref. 22) to compute bodies that support conical 
shock waves. He presents a calculation for an elliptic conical shock 
wave, and the results of the same case as computed by the method of 
Stocker. These are shown in figure 15, along with a calculation by the 
velocity-entropy program. The agreement is very good. Note that pressure 
has been plotted against the angle where \j / - arc ta^y^/z^) . 

An interesting facet of the present work is the inclusion of 
calculation of a so-called "externally iterated" density. (See eq. (43) 
and also ref. 12.) Pressures were computed at the place cp = 0 on the 
body for the case shown in figure 8 using both iterated and uniterated 
densities. The results differed by about 1.5 percent at that point. 


Comparisons With Experimental Data 


Various wind-tunnel experiments have been performed on elliptic 
conical bodies (see, e.g., refs. 23 and 24) . The shock supported by an 
elliptic conical body is not an elliptic cone, and conversely an elliptic 
conical shock wave will be supported by a body that is not elliptic in 
cross section. Thus no precise comparison can be made between results of 
tests on elliptic cones and the calculations that assume an elliptic 
conical shock wave. Some calculations were performed, however, assuming 
that the shock waves of references 23 and 24 are elliptic cones. The 
resulting bodies and pressures are shown with the experimental data in 
figures l6 and 17* Note that pressures are plotted against the angle \|r 
where V = arc ta^y^/z^) . Two attempts were made in calculating the 
conical body to compare with the actual cone tested in the = 6 
experiments of reference 24, and both results are shown in figure 17 . 


Comparison With an Approximate Theory 


A case was calculated for a circular conical shock wave with 
0 O = 30°; Moo = 10, 7 = 1.4, a = 0 , and (3 = 10 °, using the stream- function 
program. The resulting body was almost circular in cross section, and it 
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was used as the basis of a computation using the first-order yawed-cone 
theory given in reference 6 . Relevant parameters are, in the notation of 
reference 6, Mdo = 10, G s = 26 . 5 °; € = 0 . 1575 ; and- 5 = P = 0.1745* The 
surface pressure, calculated with respect to body axes, is 


p b = 0 . 217[1 + 0.584 cos (cp + 90 °)] 


(50) 


This expression is converted to shock-fixed axes by use of the relation 

tan(cp + 90°) = — + e ) - tan(cp + 90 °) (51) 

x 


where 


x - tan" 1 ^ s cos (^ + 90 D ) 


(52) 


The machine -computed surface pressure and th* approximate result are shown 
in figure l 8 . These calculations are slightly inconsistent in that the 
machine computation was carried out using 7 = 1.4, whereas the tables, 
reference 6 , use the value 7 = 1.405 • 


Systematic Variation of Parameters 


A number of cases have been calculated vith the stream-function 
program, wherein one of the initial parameters was systematically varied 
while the others were held fixed. The accuracy of these calculations was 
then checked by comparing a few cases with computations performed using 
the velocity-entropy program. Each of the t/o programs has presented 
difficulties which limit the quantity and quality of obtainable results. 
The stream- function program produces accurate results only where the 
increment A 0 -^ is small over the whole extrapolation region. This, for 
practical reasons, places lower bounds on and bi/ai and upper bounds 
on the angles a and (3. This problem has besn avoided in the velocity- 
entropy program, but the analysis of the entropy data is very time- 
consuming, and so the number of cases that can reasonably be considered 
is small. The results presented here agree closely by the two computer 
programs. 

There are many ways in which the initial parameters can be 
systematically varied. No attempt is made here to do more than sample a 
few combinations, however. The cases presented are intended mainly to 
illustrate numerical techniques. The parameter combinations that have 
been chosen are summarized in tables I to IV, along with the numbers of 
the figures where the results are plotted. 
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CONCLUDING REMARKS 


The inverse method has been applied to the problem of calculating 
bodies that support elliptic conical shock waves at finite angles of 
incidence. The entropy layer at the surface of such bodies has been 
investigated numerically, also. Two machine computing procedures were 
employed, and it is presumed that discrepancies of results are due 
primarily to differences involved in the programming of the two procedures. 

The application of two stream functions, one of which vanishes on the 
body, is a powerful technique in the solution of the three-dimensional 
conical flow problem. Future attacks on the problems of more complicated 
three-dimensional flows may employ two stream functions with advantage. 

The requirement that the shock wave be an elliptic cone makes it 
difficult to compare calculations with results of actual tests on conical 
bodies. A logical extension of the present work would be to develop 
techniques whereby the present assumed shock shapes could be perturbed 
so that a body of prescribed shape could be found in a few tries. 

The major difficulties that have been encountered in the machine 
computations could be overcome to a large extent if the segmenting 
processes of the velocity-entropy program and the two- stream-function 
formulation of the problem were incorporated into a single computing 
program. Thus good accuracy would be assured, even for very thin conical 
bodies, and the plotting of final results could be accomplished in a 
reasonably short time. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., July 26 , i960 
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TABLE I.- MACH HUMBER VARIATION AT FIXED VALUES OF bi/ai 
AND WITH 0 0 = 30° AND AT ZERO INCIDENCE 


bi/a! 

Mach no. 

Figure 

0.6 

10, 20, 100 

19 

.7 

5, 10, 20, 100 

20 

.8 

5, 10, 20, 100 

21 

.9 

3, 5, 10, 20, 100 

22 


TABLE II.- VARIATION OF 'b 1 /a. 1 AT FIXED VALUES OF 0 O AND 
WITH Moo = 10 AND AT ZERO INCIDENCE 


deg 

Wan. 

Figure 

10 

0.8, 0.9, 1.0 

23 

15 

0.7, 0.8, 0.9, 1.0 

24 

30 

0.6, 0.7, 0.8, 0.9, 1.0 

25 

40 

0.8, 0.9, l.o 

26 

50 

0.8, 0.9, l.o 

27 
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TABLE III.- ELLIPTIC CONICAL SHOCK, 0 O = 30°, AND 
WITH Moo = 10, AT FINITE INCIDENCE 


a, 

deg 

deg 

Figure 

o 

I — 1 

o o 

L T\ 

*N 

o o 

28 

29 


TABLE IV.- CIRCULAR CONICAL SHOCKS WIT I Moo = 10 AND WITH 
VARIATIONS OF 3 FOR FIXED VALUES OF 0 O 


&0) 

deg 

deg 

Figure 

10 

0 , 2.5 

30 

15 

o, 2 . 5 , 5 

31 

30 

0, 2.5, 5, 10, 20 

32 

4o 

0, 5, 10, 20 

33 

50 

o, 5, 10, 20 

34 
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START 


Load specific values of $o> 
b ( /a,, a, /3, y, Ad, and 
n for the case at hand. 


Calculate shock wave data 
using equations (25)-(28), 
(33) - (36), and (38M42). 


( ' Are any \ 

values of F( n+ 0 
negative? j 


Compute <£ -derivatives 
numerically. Calculate Pq, 
F qq, G qq using equations 

(I9)-(2I). 


Compute and print out body 
coordinates and pressure 
using equations (46-48). 


Compute /? !t and p. Print 
out p, p , F, G and S, and 
other flow-field data. 


This case is finished. Go to 
START for next case. 


Extrapolate 0 n -data to 
coordinate 9 n + \ using equa- 
tions (44) -(45). 


Figure 5-- Flow chart for the stream-function program 
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0<cf><90° for elliptic conical shock, a = /3= 0 

-90 < <f> < 90° for circular conical shock, a = 0, /3> 0 


Figure 6.- The variation of entropy with the coordinate cp in the flow 
field behind elliptic conical shocks at zero incidence, or circular 
shocks at positive angles p. 
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Surface pressure 



Figure 7*- Comparison of surface pressure in the limiting case of 

00 and 7 -+ 1 with simple Newtonian theory and the theory of 
reference 19 • 




Figure 8.- Comparison of hody shapes for an elliptic conical shock wave 
as calculated hy the two computer programs and the method of Stocker. 
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Figure 10.- The variation of the velocity conponent u in the flaw field 
behind the shock for the case where 9 0 = 30 ° > H» = 6, = 0.6, 

7 = 1.4, a = p = 0. 


A 

3 

8 

5 





Velocity, 


ko 


i 

l 


<t>, deg 


Figure 12 •- The variation of the velocity component w in the flow field 
behind the shock for the case where 9 0 ~ 30°, Moo = 6, ^x/ 3 -! = 0-6, 

7 = 1.4, a = p = 0. 





Entropy, 



Figure 13 . - The variation of entropy in the flow field behind the shock 
for the case where 0 O = 30 °, Moo = 6 , t>i/ a i = 0-6, 7 = 1*^, a = (3 - 0. 
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Figure l4.- Streamlines in the crossflow between shock and body for the 
case where 0 O = 30°, Moo = 6, ~ 0.6, 7 = 1.4, a = p = 0- 
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(b) Surface pressure 


Figure 15*- Comparison of "body shape and surface pressure for an elliptic 
conical shock wave as calculated hy the velocity-entropy program and 
the methods of Stocker and Mauger. 
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(a) Body shapes 



(b) Surface pressures 

Figure 17 •- Comparison of calculations with experiments on an elliptic 
cone of unit length having semiaxes of lengths 0.225 and 0.402. 
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(a) Body and shock wave shapes 



Figure 20 •- The effect of Mach number variation on body shape and surface 

pressure; ^>x/ a i = 0*7 series. 





Figure 21.- The effect of Mach number variation on body shape and surface 

pressure ; b -,/a-, = 0.8 series. 
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(a) Body and shock wave shapes 



(b) Surface pressure 


Figure 22.- The effect of Mach number variation on body shape and surface 

pressure; b 1 /a 1 = 0.9 series. 
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(b) Surface pressure 


Figure 23 • - The effect of the variation of on body shape and 

surface pressure; 0 O = 10° series. 
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Figure 24.- The effect of the variation of 0X1 "k°3-y shape and 

surface pressure; 0 O = 15 ° series. 
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Figure 25*- The effect of the variation of on hody shape and 

surface pressure; 0 O = 30 ° series. 
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(a) Body shape 
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(b) Surface pressure 

Figure 27-- The effect of the variation of bj/a-L on body shape and 
surface pressure; 9 0 = 50° series- 
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(b) Surface pressures 

Figure 28.- The effect of a on body shape and surface pressure; an 

elliptic conical shock wave. 
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(a) Body and shock wave shapes 
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(b) Surface pressures 


Figure 29-- The effect of £ on "body shape and surface pressure; an 

elliptic conical shock wave. 
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Figure 31.- The effect of P on "body shape and surface pressure; 
circular conical shock with 0 o = 15 ° • 



















